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An efficient numerical approach to equilibrium properties of strongly coupled sys- 
tems which include a subsystem of fermionic quantum particles and a subsystem 
of classical particles is presented. It uses an improved path integral representation 
of the many-particle density operator and allows to describe situations of strong 
coupling and strong degeneracy, where analytical theories fail. A novel numerical 
method is developed, which allows to treat degenerate systems with full account of 
the spin scatistics. Numerical results for thermodynamic properties such as inter- 
nal energy, pressure and pair correlation functions are presented over a wide range 
of degeneracy parameter. 



1 Introduction 

Thermodynamic properties of correlated Fermi systems continue to attract the 
interest of researches in many fields, including plasmas, condensed matter and 
astrophysics. All present theories have practical difficulties. Quantum kinetic 
methods, including the Kadanoff-Baym equations ErB, easily handle quantum 
and spin effects, but they become extremely complicated when correlation 
effects can not be treated pcrturbativcly. It is, therefore, important to consider 
alternative approaches, such as Monte Carlo simulations tra, as they allow for 
an efficient treatment of strong coupling phenomena. However, they still have 
difficulties in incorporating quantum and spin statistics effect. In particular, 
there has been remarkable progress in applying path integral quantum Monte 
rlo (PIMC) techniques to Bose systems H and Coulomb systems, see e.g. 
and Ref. 13 for an overview. Nevertheless, there has still been one major 
obstacle preventing efficient modelling of Fermi systems - the so-called sign 
problem resulting from summation over all permuations of the density matrix. 
In this paper, whe present a different path integral representation in which no 
sign problem appears. 
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2 Path-Integral Monte Carlo Methods 

As is well known the thermodynamic properties of a quantum system of N 
particles are fully determined by the partition function Z and, consequently, 
by the density matrix 

Z= [dqp(q,0;q,j3), p(q, 0; q',P) = (q\p\q'), (1) 



where p = exp{— f3H}, (3 — 1/kT, and q comprises the coordinates of all 
particles, q = {qi, q2, qj\r}- With an analyticaLfixpression for the density 
matrix given, one can use e.g. Monte Carlo methodsBu to evaluate the partition 
function and thermodynamic quantities. However, for a quantum system, p is, 
in general, not known, but can be constructed from its known high-temperature 
limit by means of a decomposition into n + 1 factors, each of whichjCorresponds 
to the density matrix at an n + 1 times higher temperature BHa'Ei'lj, 

p(q, 0; q', 0) = ( dq [ V . . . dq™ p(q, 0; <?« , r«) x . . . x p( q W , r< ft > ; q',f3), (2) 



where T^ i+1 ~> - = A/3 = f3/(n + 1). Eq. (|) lets one view the N-particle 
state as a loop consisting of n+1 vertices ("beads") located at intermediate 
coordinates q^ which is closed, i.e. = gW, due to the trace in Eq. ([!]). 

For quantum systems of bosons (fermions), furthermore, the spin statistics has 
to be taken into account requiring to perform an (anti-)symmetrization of the 
density matrix, after which Eq. (Q) obtains the form (for simplicity, the spin 
variables are not written explicitly) 

p(q, 0; q', (3) = -L ]T(±1) K - ( dg« . . . dq^ p(q, 0; qW , r W) 
' p v 

xp( q ^,r^;Pq',0), (3) 

where P denotes an arbitrary A^-particle permutation of the particle coordi- 
nates, and Kp is its parity. 

The unknown density matrix (|^) is efficiently computed by approximating 
the p's on the r.h.s. of Eq. (||) by the high-temperature density matrix- ,o' lT . 
For potentials which are bounded from below the simplest choice is&LfBtl 

p hT (q,T-q> , r') = p (q,r; q' , r') e^' 

p (q, T ; q i,T') = ^ w e^-^ 2 / x -, (4) 
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where po is the free-particle density matrix, U is the potential energy, and Aa 
is the thermal DeBroglie wave length corresponding to the higher temperature 
kT' = l/A/3, A A = 2tt/i 2 (t - r')/m. It is well known that for n -> oo, p hT 
converges to the exact density matrix p. Notice that the factorization of p hT 
into a kinetic and potential energy term, Eq. (^), is an approximation also. 
The error made thereby is of the order of the variation of U on the spatial 
scale A a and thus vanishes with n — > oo. This holds also for an repulsive 
Coulomb potential, whereas the attractive electron-ion interaction has to be 
represented by a bounded from below effective pair potential BtrEEl. 



3 Thermodynamic quantities of dense quantum plasmas 

The sign problem is solved by a simple transform of the intermediate electron 
coordinatesEj. To this end, we rewrite the partition function ([!]), now explicitly 
including the (classical) ion component, N e = Ni = N , 



with Q(N e ,N i) /3)= dqdrd^p(q,[r},(3). (5) 



v 

Here, the notation q is retained for the ions, [r] summarizes the electron 
coordinates with r denoting the coordinates at the beginning of the loop and 
. . . £W the dimensionless distances between neighboring vertices. Thus, 
explicitly, [r] = [r,r + A A £ (1) ,r + A A (£ (1) + £ (2) ), • . .], and q, r, £ (i jL .each are 
3iV-dimensional vectors. For the density matrix in Eq. (||) we haveB'llSl 

p(q, [r],/3) = £ Ps(q, [r], /?) = E e~^^ J[ J[ <det I^U, 

s=0 s=0 i=lp=l 

where U(q, [r], 0) = U\q) + — j- V {[/f([r], /3) + tfffo [r],/3)} , (6) 

n + 1 z — ' 

z=o 

and J7 1 denotes the ion-ion interaction energy, whereas Uf denotes the inter- 
action energy corresponding to the electronic vertex "1", and Uf l means the 
analogous electron-ion contribution. Furthermore, cj> l pp = exp[— 7r|^| 2 ] arises 
from the kinetic energy density matrix po, while 4> a b is the exchange matrix the 
elements of which are given by = exp[— 7r|(r a — rb)/\A + Ylk=i£a\ 2 ]- The 
determinant arises from the sum over the permutations, cf. Eq. (H). Finally, 
s is the number of electrons having the same spin projection. Due to the spin 
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variables, the exchange matrix contains two zero submatrices which are related 
to electrons with opposite spin projection. 

As an example of applying our method to thermodynamic properties, we 
provide the result for the equation of state, ftp = dlnQ/dV — {adlnQ /3Vda] a =i, 



(3pV 



1 



2N ' ' 6NQ 

Ni jV, 



l - £ / dq dr d£ Ps (q, [r] ,P) ]T f-r + J2 r-^ 
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\-\- \x pt \ df3® le I3e 2 {r l pt \r pt ) 

2-*>2-*> n+1 d\x pt \ + 2^2^ 
p=l t=l 1 pt ' p=lp<t 



9 M ' 



n Ni N a 



UUt^( n + VKt\ d Kt\ det|^| 



da 



(T) 



Here, $ le is the effective electron-ion pair potential, a is a scaling parameter, 



denotes the scalar product, and q p t,£, 



of two coordinate vectors: q pt = q p 
r l p - r\ + A A El=i Cpt, and 4t = r 2 



pi ) ' pt 
<7t7 r pt = r p ~ 

- qt + A A EL 



and differences 

rrt £k &k ryj 

1 1, Spi — Sp Si i ' 
1 Sp 



pt — 4"p Si i 'pi 

We underline that 



introducing the dimensionless coordinate components £ removes a substantial 
part of the complicated temperature and volume dependence from the density 
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Figure 1: Pressure of a mixture of ideal electrons and protons. Spin statistics of the electrons 
is fully included, protons are treated classically. 



matrix (Q) • This turns out to be crucial in the computation of thermodynamic 
quantities as they involve derivatives of the partition function with respect to 
V or (3. Furthermore, the anti-symmetrization propblem is reduced to comput- 
ing the determinant of (j> pt , and no explicit summation over permuations with 
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alternating sign is required. In fact, for each fixed particle configuration, the 
absolute value of each sum in curly brackets in Eq. (Q) is bounded as n — > oo, 
which is sufficient to eliminate the sign problem completely. This enables us 
to evaluate thermodynamic properties of quantum plasmas without further 
approximations on the density matrix and the integration region. For the nu- 
merical calculations, we use the standard Metropolis Monte Carlo procedure in 
which the probability of sampling configurations is proportional to \p s \, while 
the sign of p s is included into the weight function of the particle configura- 
tion. Expressions similar to (Q) are readily derived for other thermodynamic 
quantities and, along with more details, will be given elsewhere. 



4 Numerical results 

Using the above results, we have performed a series of calculations for a two- 
component electron-proton plasma. Notice that this first-principle treatment 
of the elementary plasma particles (physical picture) does fully include bound 
states. Of course, the quality of describing their properties depends on the 
potential <I )le for which we chose an effective quantum pair potential originally 
derived by Kelbg which is the exact high-temperature limit, e.g. □. Thus, 
by increasing the number of electronic vertices n, in principle, any desired 
accuracy can be achieved. To test the quality of the reproduction of quantum 
and spin statistics effects, we first consider a mixture of ideal electrons and 
protons for which the thermodynamic quantities are known analytically, e.g. 




r [a B ] 



Figure 2: Pair correlation functions for non-interacting electrons and protons at nA 3 = 4. 
The decrease of g ee for small distances reflects the Pauli principle. Also shown is the vertex 
separation distribution Pj. 
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Fig. 1 shows our numerical results for the pressure together with the theoretical 
curve. The agreement, up to values of the degeneracy parameter \ = n\ 3 
as large as 5 is remarkable. Even with only N — 16 electrons and protons 
deviations are rather small. One clearly sees that increasing the number of 
particles improves the numerical results systematically. 
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Figure 3: Pressure (upper figure) and total energy —3NkT (lower figure) of an electron- 
proton plasma for different values of the classical coupling parameter T. For comparison, 
results of other groups are showrif— iQMCNE - quantum MC without exchange, RSDWK - 
analytical model of Rieman et al. E3- 

Fig. 2 shows the correspondig e-e, i-i and e-i pair correlation functions. As 
expected, the functions g e i and gu are identical to one. The fluctuations at 
small distances reflect the maximal statistical error of the MC simulation and 
rapidly decay with increasing distance. The decay at large distances is a conse- 
quence of the periodic boundary conditions in the MC simulation. In contrast, 
the electron-electron correlation function decays at small distances reaching 
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0.5 at r — which is the expected result for particles with spin 1/2. Notice 
further the maximum of g ee which appears around the thermal wavelength A e 
and reflects weak ordering due to Fermi repulsion. To characterize the strength 
of quantum effects, we also show the probability distribution of distances be- 
tween neighboring vertices on the electron loops It has a maximum on a 
distance of the order of Aa/v^- 
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Figure 4: Pair correlation functions for a correlated electron-proton plasma with T = 0.4 
and different values of the degeneracy parameter. Notice the different length scales. Bottom 
right figure shows the vertex separation distribution P b (arbitrary units, see text). 

Let us now turn to the case of interacting electrons and protons. We have 
performed a series of calculations in which the classical coupling parameter 
r = (47rn e /3) 1 / 3 e 2 /47reoA:T was kept constant while the degeneracy was varied. 
The results for the pressure and energy are presented in Fig. 3. One can see that 
for weak coupling and small degeneracy parameters, \ < 0-5, exchange effects 
are small, and QMC simulations without exchange (open circles) are close to 
our results. However, with increasing \ an d I\ the deviations are growing 
rapidly. Even stronger are the discrepancies with analytical theories which are 
constructed as perturbations expansions, and thus are limited to small values 
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of x an d r. Most strikingly is the decrease of the energy as a function of \ 
predicted by the analytical models and QMC without exchange, which is in 
contrast to our results which show an increase for all values of T. Finally, 
Fig. 4 shows the pair correlation functions for intermediate coupling, T = 0.4 
for various degrees of degeneracy. Now, due to Coulomb repulsion, at small 
distances g ee and ga decay to zero. However, the decay of g ee is essentially 
different from that of the proton-proton function, despite the identical Coulomb 
force. The reason arc quantum exchange and tunneling effects in the electron 
subsystem which compete with the Coulomb repulsion. 
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